Comparison of markerless and marker-based motion capture systems using 95% functional limits of agreement in a linear mixed-effects modelling framework

Biomechanics analysis of human movement has been proven useful for maintenance of health, injury prevention, and rehabilitation in both sports and clinical populations. A marker-based motion capture system is considered the gold standard method of measurement for three dimensional kinematics measurements. However, the application of markers to anatomical bony points is a time consuming process and constrained by inter-, intra-tester and session reliability issues. The emergence of novel markerless motion capture systems without the use of reflective markers is a rapidly growing field in motion analysis. However an assessment of the level of agreement of a markerless system with an established gold standard marker-based system is needed to ensure the applicability of a markerless system. An extra layer of complexity is involved as the kinematics measurements are functional responses. In this paper a new approach is proposed to generate 95% functional limits of agreement (fLoA) using the linear mixed-effects modelling framework for hierarchical study designs. This approach is attractive as it will allow practitioners to extend their use of linear mixed models to assess agreement in method comparison studies in all domains where functional responses are recorded.

is the variance function assumed to be continuous over δ.The choice of the function g(.) depends on the context.If it is believed that the variance of the within-group error increases linearly with time then the variance model would be and the corresponding variance function is There are different variance functions and correlation structures available in the nlme package in R (R Core Team, 2021; Pinheiro and Bates, 2000).This makes variance functions and correlation structures easily accessible when fitting a mixed-effects model in R. Table 1 lists the different variance functions available in the nlme package in R.
In the mixed-effects modelling framework, the correlation structure is used to model the dependency of the within-group error.Historically, these dependency structures have been developed in two areas of statistics: time series analysis and spatial data analysis.The difference between these two types of data is that response variable in time series is usually indexed by one variable (i.e.time) whereas in spatial data, the response is usually indexed by two coordinates of a spatial plane.As the motion capture study involves a time series type response, dependency structures needed for time series analyses will be considered.
To develop a general structure similar to the variance function, consider the case where the dependency of the within-group errors only depend on some position vector p ij .In the situation where a univariate response is considered, this position vector is just a scalar.In the case of spatial data analysis for example, this position vector may contain a multidimensional vector.It is also assumed that the correlation structure is isotropic.This means the correlation between two errors depends only through some distance, say d(p ij , p ik ), between these two positional vectors relating to those errors.
A general correlation structure for errors is therefore as follows (Pinheiro and Bates, 2000): where, ρ is the vector of correlation parameters and h(.) is the correlation function.
This is a very general structure for modelling the dependency in the within-group errors both for the spatial data and time series data.Since the nature of the data in this thesis is similar to time series data, from now on the focus will be on correlation structures relevant for time series data.The correlation in time series data is known as serial correlation.Since only time series data are considered, p ij will only contain a scalar position index which will be denoted as p ij since it is no longer a vector.
The isotropic assumption will further be simplified to a situation where the correlation only depends on the absolute value of the difference between two position indexes.The general serial correlation structure can now be modelled as (Pinheiro and Bates, 2000) cor In time series data, the correlation function h(.) is referred to as autocorrelation function.A nonparametric estimate of the autocorrelation function is known as the empirical autocorrelation function and can be used to examine the serial correlation in the data.Let r ij = (y ij − ŷij )/σ ij denote the standardised residual from a fitted mixed-effect model where σ2 ij is the estimate of Var(ϵ ij ) = σ 2 ij , then the empirical autocorrelation function at lag l is defined as (Pinheiro and Bates, 2000) ρ where, n is the total number of subjects, n i is the number of observations for the i th subject, N (l) is the number of residual pairs used pairs used in the summation to define the numerator of ρ(l), and N (0) is the total number of residuals.
The simplest serial correlation structure is compound symmetry, which can be defined as follows: In this situation the autocorrelation function is h(l, ρ) = ρ, where l = |j − k|.This is a very simplistic correlation structure and might not be very useful in general.The other extreme is the general correlation structure, where the autocorrelation function is This may also not be a useful correlation structure as it requires many correlation parameters to be estimated.Therefore, this may be useful only to find a more parsimonious correlation structure for exploratory purposes.
The correlation structure that is most relevant for this thesis comes from a different class of linear stationary models: autoregressive models and moving average models.These models assume that the measurements were taken at discrete time points.Consider ϵ t as the measurement at time point t.The distance, or lag, between two measurements ϵ t and ϵ s is |t − s| where lag-1 means the measurements are one unit apart.The autoregressive model assumes that the measurement at the current time is linearly dependent upon the previous measurements plus homoscedastic white noise, a t , centred at zero, E(a t ) = 0 (Pinheiro and Bates, 2000).
The number of previous observations on which the current measurement depend upon is called the order of the autoregressive model.The order of the autoregressive model here is p and the model is denoted as an AR(p) model.Note that p is used previously to denote the number of fixed-effects in a LMM.For this section, p will be used as the order of an AR model.The model also contains the same number of correlation parameters, ϕ = (ϕ 1 , . . ., ϕ p ) ′ as the order of the model.The correlation function for an AR(1) model is follows (Pinheiro and Bates, 2000): where k is the distance between two time points.The correlation function beyond the AR(1) model does not have any simple representation.It is defined recursively through the difference equation (Pinheiro and Bates, 2000) Moving average models assume that the current observations are a linear combination of independent and identically distributed white noise terms (Pinheiro and Bates, 2000) The number of white noise terms with lag, q, is the order of the moving average model which is denoted by MA(q) model.Note that q is used previously to denote the number of random-effects in a LMM.
For this section, q will be used as the order of a MA model.There are q correlation parameters in this model θ = (θ 1 , . . ., θ q ) ′ .
The correlation function for an MA(q) model for the observations with k distance apart is as follows A combination of an autoregressive and a moving average model is called an autoregressive-moving average model and denoted by ARMA (p, q) with the order p for the autoregressive model and order q for the moving average model.This model can be written as follows (Pinheiro and Bates, 2000): The correlation function for this model is defined recursively as follows (Pinheiro and Bates, 2000): ϕ 1 h(|k − 1|, ρ) + . . .+ ϕ p h(|k − p|, ρ)+ θ 1 ψ(k − 1, ρ) + . . .+ θ q ψ(k − q, ρ), k = 1, 2, . . ., q ϕ 1 h(|k − 1|, ρ) + . . .+ ϕ p h(|k − p|, ρ), k = q + 1, q + 2, . . ., where ψ(k, ρ) = E(ϵ t−k a t )/Var(a t ).Table 2 lists the available correlation functions in the nlme package.Note that this table includes correlation structures for time series data and correlation structures for spatial data.Only the correlation structure for time series data has been discussed in this thesis and the details for the spatial correlation structure can be found in Pinheiro and Bates (2000).

Table 2 :
Name of the correlation functions available in the nlme package in R.